# Analysis the relationship between TF expression and efficacy
# 2024/7/17

# + Init ----
library(readxl)
library(tidyverse)
library(dplyr)
library(ggsci)
library(ggplot2)
library(ggpubr)
library(forcats)
library(ggblanket)
library(ggrepel)
library(ggtrendline)
# library(export)

rm(list = ls())

tumor2en <- function(tumor.cn) {  
  tumor.en <- tumor.cn
  
  # CC:宫颈癌, PC:胰腺癌, UC:尿路上皮癌, OC:卵巢癌, HNSCC:头颈部鳞癌, 
  #   FTC:输卵管癌, Prostate:前列腺癌
  # 为了分析方便，此处输卵管癌也算作卵巢癌
  tumor.en <- str_replace(tumor.en, "宫颈癌", "CC")
  tumor.en <- str_replace(tumor.en, "胰腺癌", "PC") # Pancreatic cancer
  tumor.en <- str_replace(tumor.en, "尿路上皮癌", "UC")
  tumor.en <- str_replace(tumor.en, "卵巢癌", "OC")
  tumor.en <- str_replace(tumor.en, "头颈鳞癌", "HNSCC")
  tumor.en <- str_replace(tumor.en, "鼻咽癌", "NPC") # nasopharyngeal carcinoma
  # tumor.en <- str_replace(tumor.en, "输卵管癌", "FTC") # Fallopian tube cancer
  tumor.en <- str_replace(tumor.en, "输卵管癌", "OC") # Fallopian tube cancer
  tumor.en <- str_replace(tumor.en, "前列腺癌", "Prosate")
  
  return(tumor.en)
}


# + Load data -----
# state2en <- function(state.cn) {
#   state.en <- state.cn
#   state.en <- str_replace(state.en, "终止", "Discontinue")
#   state.en <- str_replace(state.en, "治疗", "Ongoing")

#   return(state.en)
# }

# 为了按疗效排序
resp.level <- c("CR", "PR", "SD", "PD", "NE", "NA")

# Load patient info
data.file = "D:/Work/TF/Ph1/X28_explore.xlsx"

pts.data = read_excel(data.file, sheet = "EffectPredictor")
pts.data$Tumor <- tumor2en(pts.data$Tumor) # Translate into EN

# Transfer to factor
pts.data <- pts.data |> 
  select(Patient, CurrentState, Sex, Age, Tumor, Dose, BOR, NLR_mean, mH_score, mTF1, mTF2, mTF3) |> 
  mutate(CurrentState = factor(CurrentState, levels = c("治疗", "终止"), 
                               labels = c("Ongoing", "Discontinue"))) |> 
  mutate(Sex = factor(Sex, levels = c("男", "女"), labels = c("M", "F"))) |> 
  mutate(BOR = factor(BOR, levels = resp.level, ordered = TRUE))

# Drop patient data without TF and BOR
tf.data <- pts.data |> 
  filter(!is.na(mH_score)) |> 
  filter(mH_score != "NA") %>% 
  mutate(BOR = factor(BOR, levels = resp.level, ordered = TRUE)) |> 
  filter(BOR != "NA") |> 
  select(Patient, Sex, Age, Tumor, Dose, BOR, NLR_mean, mH_score, mTF1, mTF2, mTF3)

# tf.data <- tf.data %>% 
#   mutate(H_Resp = )

tf.data$mH_score <- as.numeric(tf.data$mH_score)

glimpse(tf.data)


# + TF expression -----

## 1. Tumor and TF
# ggpubr
ggboxplot(tf.data, x = "Tumor", y = "mH_score", 
          label = "mH_score",
  palette = "jco", color = "Tumor",
  # add = "dotplot",
  xlab = "", ylab = "H score"
)

# ggplot2 效果同上，但根据TF表达量降序排列
tumor.boxplot <- tf.data |>
  mutate(Tumor = fct_reorder(Tumor, desc(mH_score), .fun = "median")) |>
  ggplot(aes(x = Tumor, y = mH_score, fill = Tumor)) +
  geom_boxplot(alpha = 0.5) +
  # geom_label(tf.data, aes(label = mH_score)) +
  geom_point(aes(color = BOR), size = 4, alpha = 0.7) + # geom_jitter()
  # geom_jitter(size = 2, alpha = 0.7) +
  labs(x = "", y = "H score") +
  theme_bw() +
  theme(legend.position = 'none') +
  scale_fill_lancet()

tumor.boxplot +
  # geom_text_repel(aes(label = mH_score), size=3) +
  geom_text_repel(aes(label = BOR), size = 3, alpha = 0.8)

ggsave("tf_resp.png", width = 8, height = 6, dpi = 600)

# # ggblanket
# gg_boxplot(data = tf.data, x = Tumor, y = mH_score,
#            col = Tumor)

## 2. BOR and TF
tf.data |> 
  ggboxplot(x = "BOR", y = "mH_score", 
    palette = "jco", color = "BOR", add = "dotplot",
  title = "All cases", xlab="", ylab="H score")

# ggplot2
tf.data %>% 
  ggplot(aes(x = BOR, y = mH_score, fill = BOR)) +
  geom_boxplot(alpha=0.7) +
  geom_point(size = 2) +
  geom_text_repel(aes(label = Tumor), alpha = 0.6) +
  labs(x = "", y = "H score") +
  theme_bw() +
  theme(legend.position = 'none') + # without Legend
  scale_fill_jama()

ggsave("tf_bor.png", width = 8, height = 7, dpi = 600)

# Sub group
tf.data |> 
  filter(Dose != 0.6) |> 
  ggboxplot(x = "BOR", y = "mH_score", 
    palette = "jco", color = "BOR", add = "jitter",
  title = "Without 0.6 mg/kg")

tf.data |> 
  filter(mTF3 > 0) |> 
    ggboxplot(x = "BOR", y = "mH_score", 
  palette = "jco", color = "BOR", add = "jitter",
  title = "TF 3+", xlab = "")

  tf.data |> 
    filter(mTF3 + mTF2 > 0) |> 
    ggboxplot(x = "BOR", y = "mH_score", 
      palette = "jco", color = "BOR", add = "jitter",
      title = "TF 2+", xlab = "")

#+ TF Linear regression ----
tumor.tf <- tf.data |>
  group_by(Tumor) |>
  summarise(mH_score = mean(mH_score))

# Target lesion data
resp.file = "D:/Work/TF/Ph1/XNW28012_SubjectList.xlsx"
recist.df <- read_xlsx(resp.file,
  sheet = "TargetLesion",
  range = cell_cols("A:E"),
  col_types = c("text", "date", "numeric", "numeric", "text"))

recist.df <- recist.df %>% 
  filter(Sum > 0) # Del NA (Del patients without target lesions)
# Baseline
baseline.df <- recist.df %>% 
  filter(No.eva == 0)

# Get the sum of target lesions of the best response
best.resp <- recist.df %>%
  filter(No.eva > 0) %>%
  group_by(Patient) %>%
  summarize(Min = min(Sum))

# Calculate the Max change from baseline
#   then join the patients info
shrinkage.df <- left_join(best.resp, baseline.df, by = "Patient") %>%
  mutate(MaxChange = - (Min - Sum) / Sum * 100) %>%
  filter(!is.na(MaxChange)) |> 
  select(Patient, MaxChange)

change.tf <- left_join(tf.data, shrinkage.df, by = "Patient")
change.tf <- na.omit(change.tf)

# Regression
fit <- lm(MaxChange ~ mH_score, data = change.tf)

summary(fit)

plot(change.tf$mH_score, change.tf$MaxChange, xlab = "H score", ylab = "Max Shrinkage", main = "Max shrinkage")
abline(fit)

# Regression plot
ggplot(change.tf, aes(x = mH_score, y = MaxChange)) +
  geom_point() +
  geom_smooth(method = "lm") +
  theme_bw() +
  scale_color_bmj()

# Better than above (ggplot2)
library(ggtrendline)
ggtrendline(change.tf$mH_score, change.tf$MaxChange, 
            CI.fill = NA, model = "line2P",
            linewidth = 1.0,
            xlab = "H score", ylab = "Max Shrinkage%") +
  geom_point(aes(change.tf$mH_score, change.tf$MaxChange),
             color = "steelblue", shape = 1, size = 2) +
  theme_bw()

ggsave("tf_regr.png", width = 8, height = 6, dpi = 600)

# fit2 <- lm(MaxChange ~ mTF1+mTF2+mTF3, data = change.tf)
# summary(fit2)
# plot(change.tf$mTF1, change.tf$MaxChange, xlab = "H score", ylab = "Max Shrinkage", main = "Max shrinkage")
# abline(fit2)


# + NLR ----
nlr.boxplot <- tf.data %>% 
  ggplot(aes(x = BOR, y = NLR_mean, fill = BOR)) +
  geom_boxplot(alpha = 0.7) +
  geom_point() +
  geom_text_repel(aes(label = Patient), size = 3) +
  scale_fill_nejm() +
  theme_bw() +
  theme(legend.position = "none") +
  labs(x="", y="NLR")

nlr.boxplot

# ggsave("nlr.png", width = 8, height = 6, dpi = "retina")

fit <- lm(MaxChange ~ NLR_mean, data = change.tf)
summary(fit)

plot(change.tf$NLR_mean, change.tf$MaxChange, 
     xlab = "H score", ylab = "Max Shrinkage", main = "Max shrinkage")
abline(fit)

ggtrendline(change.tf$NLR_mean, change.tf$MaxChange, 
            CI.fill = NA, model = "line2P",
            linewidth = 1.0,
            xlab = "NLR", ylab = "Max Shrinkage%") +
  geom_point(aes(change.tf$NLR_mean, change.tf$MaxChange),
             color = "steelblue", shape = 1, size = 2) +
  theme_bw()
